# -*- coding = utf-8 -*-
import pandas as pd
import numpy as np
import math

if __name__ == "__main__":

    # 输入要求解的矩阵
    matrix = pd.DataFrame(
        np.array([
            [-3, -2, 0, 0, 0],
            [1, 1, 1, 0, 40],
            [2, 1, 0, 1, 60]
        ]),
        index=['obj', 'x2', 'x3'],
        columns=['x0', 'x1', 'x2', 'x3', 'b']
    )
    # 下标集 表示x0,x1,x2,x3
    subscribe_set = [i for i in range(matrix.shape[1] - 1)]
    # subscribe_set: [0, 1, 2, 3]
    A = matrix.iloc[1:, :-1]
    m = A.shape[0] - 1  # 方程个数
    n = A.shape[1] - 1  # 未知量个数

    # 1.选择初始解
    x = np.array([0, 0, 40, 60]).reshape((4, 1))
    B = np.array([[1, 0],
                  [0, 1]])
    IB = np.array([2, 3])  # 基变量索引
    NB = np.array([0, 1])  # 非基变量索引
    # ------------------*--------------------
    # 2.计算单纯形乘子 p
    CB = np.zeros((2, 1))
    for i in range(m + 1):
        CB[i] = matrix.iloc[0, IB[i]]
    p = np.zeros((2, 1))
    p = np.dot(np.linalg.inv(B.T), CB)
    p_T = p.reshape((1, 2))

    # 3.计算即约费用
    r = np.zeros((1, 2))
    for i in range(m + 1):
        A_i = np.array(matrix.iloc[1:, NB[i]])
        A_i = A_i.reshape((2, 1))
        c_i = np.array(matrix.iloc[0, NB[i]])
        r_i = c_i - np.dot(p_T, A_i)
        r[0, i] = r_i
        # r= [[-3. -2.]]

    # 4.检查最优性
    r_list = []
    for i in range(2):
        r_list.append(int(r[0, i]))
    # r_list = [-3, -2]
    flag = 0
    for i in r_list:
        if i < 0:
            flag += 1
        else:
            continue
    # ---------------------------------------
    # ------------------2--------------------
    itera = 0
    while flag != 0:
        itera += 1
        print(itera)
        # 2.计算单纯形乘子 p
        CB = np.zeros((2, 1))
        for i in range(m + 1):
            CB[i] = matrix.iloc[0, IB[i]]
        p = np.zeros((2, 1))
        p = np.dot(np.linalg.inv(B.T), CB)
        p_T = p.reshape((1, 2))

        # 3.计算即约费用
        r = np.zeros((1, 2))
        for i in range(m + 1):
            A_i = np.array(matrix.iloc[1:, NB[i]])
            A_i = A_i.reshape((2, 1))
            c_i = np.array(matrix.iloc[0, NB[i]])
            r_i = c_i - np.dot(p_T, A_i)
            r[0, i] = r_i
            # r= [[-3. -2.]]

        # 4.检查最优性
        r_list = []
        for i in range(2):
            r_list.append(int(r[0, i]))
        # r_list = [-3, -2]
        flag = 0
        for i in r_list:
            if i < 0:
                flag += 1
            else:
                continue
        if flag == 0:
            print('找到最优解')
            # 更新IB
            IB_b = list(IB)
            IB = np.delete(IB, IB_b.index(out_idx))
            IB = np.append(IB, q)
            print('最优解为：')
            for i in IB:
                print("i=",i)
                print('x', i + 1, ':', x[i])
            print('obj=', -3 * x[0] + (-2) * x[1])
        if flag != 0:
            print('现行解不是最优解，选择入基变量：')
            # 5.选择入基变量
            q = in_idx = r_list.index(min(r_list))  # 入基变量索引
            in_x = NB[in_idx]  # 入基变量
            # in_x 0
            # in_idx 0

            # 6.转移方向
            d = np.zeros((2, 1))
            A_i = matrix.iloc[1:, NB[in_idx]]
            d = -np.dot(np.linalg.inv(B), A_i)
            # d=[-1. -2.]
            d = d.reshape((2, 1))
            # 转移方向为d_q
            d_q = np.zeros((1, n + 1))
            d_q = d_q.reshape((4, 1))
            d_q[q, 0] = 1
            d_q[IB[0]] = d[0]
            d_q[IB[1]] = d[1]
            for NB_not_inB in NB:
                if NB_not_inB != q:
                    d_q[NB_not_inB] = 0

            flag2 = 0
            for i in d_q:
                if i >= 0:
                    flag2 += 1
                else:
                    continue
            if flag2 == 4:
                print('解无界')
            else:
                print('没有检查出无界性，继续')

            # 7.转轴
            lamda = 0
            lamda_list = []
            lamda_list_0 = []
            for i in range(n + 1):
                if d_q[i] < 0:
                    x_i = x[i]
                    d_i = d_q[i]
                    lamda_list.append(-x_i / d_i)
                    lamda_list_0.append(-x_i / d_i)
                else:
                    lamda_list_0.append(0)
            lamda = min(lamda_list)

            # 选择出基
            for i in range(n + 1):
                if lamda_list_0[i] == lamda:
                    out_idx = i
                    break

            # 8.更新
            # 更新x
            for i in range(n + 1):
                if i == in_idx:
                    x[i] = lamda
                if i in IB:
                    x[i] = x[i] + lamda * d_q[i]
                else:
                    continue
            # 更新B
            A_in = np.array(matrix.iloc[1:, q])
            A_in = A_in.reshape((2, 1))
            A_out = np.array(matrix.iloc[1:, out_idx])
            A_out = A_out.reshape((2, 1))
            e_p = np.zeros((1, 2))
            for i in range(2):
                if IB[i] == out_idx:
                    e_p[0, i] = 1
            e_p = e_p.reshape((1, 2))
            B = B + np.dot((A_in - A_out), e_p)
            # 更新IB
            print("in:", q)
            print("out:", out_idx)
            # IB
            IB_b = list(IB)
            IB = np.delete(IB, IB_b.index(out_idx))
            IB = np.append(IB, q)

            # 更新NB
            NB_b = list(NB)
            NB = np.delete(NB, NB_b.index(q))
            NB = np.append(NB, out_idx)
